Demographics and reproductive biology of H. schistosus from bycatch on the west coast of India

Analysis for Manuscript

Demographics of H. schistosus across time

Sample size

The number of snakes (N) sampled from bycatch of trawlers, shore seines and gillnetss.

snakes%>%
  group_by(Species)%>%
  count(name = "N")
Species N
Hydrophis curtus 282
Hydrophis schistosus 967

The number of trips (n) and fishing effort (haul-hours) sampled for sea snakes in the current study.

snakes%>%
  filter(Gear.Type != "")%>%
  select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
  distinct()%>%
  # Calculating tow duration
  mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
                                        is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
  group_by(Gear.Type)%>%
  summarise(n = n(), # counting number of trips sampled
            haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort
Gear.Type n haul.hours
GillNet 340 535.9633
Rampan 46 190.4500
Trawler 76 285.1600

Age structure

# calculating the mean SVL across years

age.yr <- snakes%>%
  group_by(Species, Year)%>%
  summarise(N = n(),
            mean = mean(Snout.to.Vent..cm., na.rm = T))

# SVL cut off for age classes

maturtity <- snakes%>%
  group_by(Species)%>%
  count()%>%
  mutate(juv = 35,
         adult = ifelse(Species == "Hydrophis curtus", 54, 65))

Distribution fo SVL across years

# plotting distribution of SVL across years

snakes%>%
  filter(Snout.to.Vent..cm. > 20)%>%
  ggplot(aes(Year, Snout.to.Vent..cm.))+
  geom_violin(fill = "grey")+
  geom_boxplot(width = 0.1)+
  geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
  geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
  #stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
  geom_label(data = age.yr, aes(Year, 10, label = N))+
  facet_wrap(~Species, ncol = 1, scales = "free_y")+
  labs(y = "Snout to vent length (cm)")+
  theme(strip.text = element_text(face = "italic"))

ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)

We did not find any H. curtus neonates in our sampling. Majority of H. curtus were juveniles. H. schistosus populations consisted of mostly adults.

Proportion of developmental classes across years

pop.yr <- snakes%>%
  filter(!is.na(Class))%>% # fix 2016 svl
  group_by(Species, Year)%>%
  # Counting number of individuals in each developmental class
  count(Class)%>%
  complete(Class, fill = list(n = 0))%>%
  # Total number sampled in each year
  mutate(N = sum(n))%>%
  group_by(Year, Species, Class)%>%
  # Proportion of each developmental class
  mutate(prop.age = n/N)

# Summary stats

pop.yr%>%
  group_by(Species, Class)%>%
  skimr::skim(prop.age)%>%
  skimr::yank("numeric")%>%
  select(c(Species, Class, mean, sd))

Variable type: numeric

Species Class mean sd
Hydrophis curtus Adult 0.43 0.31
Hydrophis curtus Juvenile 0.57 0.31
Hydrophis curtus Neonate 0.02 NA
Hydrophis schistosus Adult 0.88 0.09
Hydrophis schistosus Juvenile 0.14 0.07
Hydrophis schistosus Neonate 0.03 0.02

Testing change in proportion of developmental classes across years

pop.yr%>%
  filter(Year != 2019)%>%
  group_by(Species, Class)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for proportion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(Class:p.value)
Species Class estimate1 estimate2 estimate3 statistic p.value
Hydrophis curtus Adult 0.8750000 0.3125000 0.3551402 9.2546983 0.0097807
Hydrophis curtus Juvenile 0.1250000 0.6875000 0.6261682 8.8509715 0.0119684
Hydrophis curtus Neonate NA NA NA 97.2336449 0.0000000
Hydrophis schistosus Adult 0.8846154 0.7894737 0.8613445 5.9949721 0.0499124
Hydrophis schistosus Juvenile 0.0769231 0.2105263 0.1239496 9.7219255 0.0077430
Hydrophis schistosus Neonate 0.0384615 0.0147059 NA 0.0189777 0.8904305

Age structure does not change significantly over a four year period from 2016 to 2019.

Change in SVL distribution over years

library(car)

snakes%>%
  filter(Year != "2019")%>%
  group_by(Species)%>%
  select(Year, Snout.to.Vent..cm.)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
          sumr = map(mod, broom::tidy),
          stat = map(mod, car::Anova))%>%
  select(stat)%>%
  unnest()
Species Sum Sq Df F value Pr(>F)
Hydrophis schistosus 9866.294 2 21.86203 0.00e+00
Hydrophis schistosus 168785.548 748 NA NA
Hydrophis curtus 2641.334 2 10.07170 7.05e-05
Hydrophis curtus 24258.408 185 NA NA

There was a significant change in SVL distribution of species across years.

month.svl <- HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month)%>%
  group_by(Month)%>%
  # Calculating mean SVL in each month
  summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))

# Month of observed births

births <- data.frame(Species = "Hydrophis schistosus", Month = "April")
# plotting distribution of SVL across months

HS%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(nesting(Species), Month)%>%
  droplevels()%>%
  ggplot(aes(Month, Snout.to.Vent..cm.))+
  geom_violin(fill = "light grey")+
  #geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
  geom_boxplot(width = 0.1)+
  geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Snout to vent length (cm)")

Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.

Proportion of jeuveniles in each month of 2018

snakes%>%
  filter(Year == "2018", # data not sufficient for other years
         Snout.to.Vent..cm. > 20)%>% # removing erroneous data
  mutate(Month = factor(Month, levels = month.name))%>%
  group_by(Species, Month)%>%
  summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
  spread(Month, prop.neonate)
Species January February March April May October November December
Hydrophis curtus 0 0 0.50 NA NA 0 0.0769231 0
Hydrophis schistosus 0 0 0.02 0.05 0.0555556 0 0.0000000 0

Sex ratios

Plotting sex ratios for each species across years

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by( Species, Year)%>%
  summarise(N = n(),
            females = sum(Sex == "Female"),
            prop.female = females/N)%>%
  ggplot(aes(Year, prop.female, fill = Species))+
  geom_col(width = 0.5, col = "black", position = position_dodge())

Proportion of females to entire population stays constant over the sampling period.

Proportion of females encountered by species

snakes%>%
  group_by(Species)%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())
Species females N
Hydrophis curtus 102 212
Hydrophis schistosus 304 618

Testing shift in sex ratios of sea snakes in bycatch

# is proportion of females different from 0.5?

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>%
  group_by(Species)%>%
  summarise(females = sum(Sex == "Female"),
            N = n())%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)), # Z - test for proportion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(-c(method, alternative))
Species estimate statistic p.value parameter conf.low conf.high
Hydrophis curtus 0.4811321 0.2311321 0.6306857 1 0.4125066 0.5504525
Hydrophis schistosus 0.4919094 0.1310680 0.7173273 1 0.4518628 0.5320580

Sex ratio is not different from 0.5; p = 0.71.

Mortality in bycatch

Overall mortality rate

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species n N prop.dead
Hydrophis curtus 120 282 0.4255319
Hydrophis schistosus 174 967 0.1799380

Testing difference in mortality across species

snakes%>%
  group_by(Species)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  ungroup()%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for propotion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)%>%
  rename(`H. schistosus` = estimate2,
         `H. curtus` = estimate1)
H. curtus H. schistosus statistic p.value
0.4255319 0.179938 71.81006 0

H. curtus had as significatly higher mortality rate in bycatch than H. schistosus

Mortality rate by age class

Are juveniles more vulnerable to bycatch mortality than adults?

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  #calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Class n N prop.dead
Hydrophis curtus Adult 40 62 0.6451613
Hydrophis curtus Juvenile 41 125 0.3280000
Hydrophis schistosus Adult 138 648 0.2129630
Hydrophis schistosus Juvenile 5 105 0.0476190
Hydrophis schistosus Neonate 3 8 0.3750000

Testing difference in mortality rates across age classes within species

snakes%>%
  filter(!is.na(Class))%>%
  group_by(Species, Class)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Caluclating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)), # Z - test for proportion
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
Species estimate1 estimate2 statistic p.value
Hydrophis curtus 0.6451613 0.328000 15.71186 0.0000738
Hydrophis schistosus 0.2129630 0.047619 17.68174 0.0001447

H. curtus had the highest mortality rate overall. H. schistosus juveniles had the lowest.

Mortality rate by sex

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Species Sex n N prop.dead
Hydrophis curtus Female 36 102 0.3529412
Hydrophis curtus Male 52 110 0.4727273
Hydrophis schistosus Female 62 304 0.2039474
Hydrophis schistosus Male 72 314 0.2292994

Testing difference in mortality rates across sexes withing species

snakes%>%
  filter(Sex != "")%>%
  group_by(Species, Sex)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
Species estimate1 estimate2 statistic p.value
Hydrophis curtus 0.3529412 0.4727273 2.6538718 0.1032980
Hydrophis schistosus 0.2039474 0.2292994 0.4448479 0.5047918

There was no difference in mortality rates across sexes in either species.

Mortality and female reproductive status

Are gravid females more susceptible to bycatch mortality than the rest of our sample?

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  # Calculating mortality rate
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)
Gravid n N prop.dead
41 151 0.2715232
Yes 16 85 0.1882353

Testing thing difference in mortality rate of H. schistosus females by reproductive state

HS%>%
  filter(Sex == "Female",
         Class == "Adult")%>%
  group_by(Gravid)%>%
  count(Condition.at.encounter..D.A.)%>%
  mutate(N = sum(n),
         prop.dead = n/N)%>%
  filter(Condition.at.encounter..D.A. == "Dead")%>%
  select(-Condition.at.encounter..D.A.)%>%
  ungroup()%>%
  nest()%>%
  mutate(test = map(data, ~prop.test(.$n, .$N)),
         sumr = map(test, broom::tidy))%>%
  select(sumr)%>%
  unnest()%>%
  select(estimate1:p.value)
estimate1 estimate2 statistic p.value
0.2715232 0.1882353 1.629856 0.2017229

Observed reproductive cycle of H. schistosus

Proportion of gravid females in the sample

# percentage of gravid females in sample

HS%>%
  count(Gravid)%>%
  mutate(N = sum(n))%>%
  mutate(prop.gravid = n/N)%>%
  filter(Gravid == "Yes")%>%
  select(-c(Gravid))
n N prop.gravid
88 967 0.0910031

Number of gravid females encountered per year

# checking the number of gravid females per year

HS%>%
  group_by(Year)%>%
  filter(Gravid == "Yes")%>%
  count(Gravid)%>%
  spread(Gravid, n)%>%
  rename(`n(Gravid Females)` = Yes)
Year n(Gravid Females)
2016 1
2018 75
2019 12

Monthly variation in proportion of gravid females

Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.

# calculating the proportion of gravid females per month

gravid <- HS%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
  group_by(Month)%>%
  summarise(N = n(),
            gravid = sum(Gravid == "Yes"),
            prop.gravid = gravid/N)

gravid
Month N gravid prop.gravid
January 102 12 0.1176471
February 134 34 0.2537313
March 129 24 0.1860465
April 54 7 0.1296296
May 36 1 0.0277778
October 16 0 0.0000000
November 69 3 0.0434783
December 34 6 0.1764706

Describing change in the relative proportions of gravid females across months and years of sampling.

# plotting prop gravid per month

gravid%>%
  mutate(Month = factor(Month, levels = month.name))%>%
  complete(Month, fill = list(prop.gravid = 0))%>%
  ggplot(aes(Month, prop.gravid))+
  geom_point(size = 3)+
  geom_path(aes(group = 1), size = 1, linetype = "dashed")+
  geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
               arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
  geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
  geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
  geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
  geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
  scale_x_discrete(guide = guide_axis(n.dodge = 2))+
  labs(y = "Proportion of gravid females")

Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.

Development of embryos and eggs

Sample size

embryos%>%
  summarise(n.mothers = length(unique(Field.Code)),
            n.embryos = length(unique(Embryo.Code)))
n.mothers n.embryos
29 235

Summary statistics of egg measurements

embryos%>%
  select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
  skimr::skim()%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable mean sd
Egg.Length..mm.. 39.71 11.72
Egg.Width..mm.. 32.23 9.01
Egg.Weigth..g.. 14.30 5.45
Snout.to.Vent..cm. 17.72 5.24
Embryo.Weight..g. 4.83 6.28

Sex ratios within clutches

embryos%>%
  group_by(Field.Code)%>%
  filter(Sex != "")%>%
  summarise(prop.female = sum(Sex == "Female")/n())%>%
  skimr::skim(prop.female)%>%
  skimr::yank("numeric")%>%
  select(c(skim_variable, mean, sd))

Variable type: numeric

skim_variable mean sd
prop.female 0.53 0.31

Number of feamle embryos

embryos%>%
  filter(Sex != "")%>%
  summarise(females = sum(Sex == "Female"),
            N = n())
females N
86 166

Testing clutch sex ratio

broom::tidy(prop.test(86, 166, p = 0.5))%>% # Z - test
  select(estimate:conf.high)
estimate statistic p.value parameter conf.low conf.high
0.5180723 0.1506024 0.6979603 1 0.4395567 0.5957383

The sex ratio in clutches is not significantly different from 0.5.

Infertility rate

Percentage of eggs with no embryos present:

embryos%>%
  count(Embryo)%>%
  mutate(N = sum(n), rate = n*100/N)%>%
  select(-N)
Embryo n rate
Absent 5 2.12766
Present 230 97.87234

Matrotrophy

Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?

Change in yolk weight with embryo weight

embryos%>%
  ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
  geom_point(size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  labs(x = "Embryo weight (g)", y = "Yolk weight (g)")

Linear model to test change in yolk weight with embryo development:

embryos%>%
  select(Yolk.Weight..g., Embryo.Weight..g.)%>%
  nest()%>%
  # Yolk weight is normally distributed
  mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.0905565 0.2965229 44.14687 0 0.4071222
Embryo.Weight..g. -0.4526283 0.0450509 -10.04705 0 0.4071222

Yolk weight decreases by 0.45 g for every 1g increase in embryo weight.

Change in egg mass with embryo development

embryos%>%
  ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
  geom_point(aes(col = Yolk.Weight..g.), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_color_viridis(name = "Yolk Weight (g)")+
  labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")

Linear model to test relation ship between total egg mass and embryo development:

embryos%>%
  select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
  nest()%>%
  # Total egg mass is normally distributed
  mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 13.5441434 0.4114909 32.91481 0 0.2188866
Embryo.Weight..g. 0.3676922 0.0557912 6.59050 0 0.2188866

Total egg mass increases by 0.36 g for every 1g increase in embryo weight.

Note: Embryo weight is log-normally distributed.

Along with decrease of fat bodies of gravid females (Figure A1) as embryos develop indicates matrotrophic nutrition.

Reproductive effort

Does the amount of enery expended by female H. schistosus reduce with age?

# Cleaning reproductive effort data

re <- embryos%>%
  left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
  select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
         Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
         -Body.Length..cm..x, - Body.Length..cm..y)%>%
  filter(Species == "Hydrophis schistosus")%>%
  rename(Embryo.SVL = Snout.to.Vent..cm..x,
         Embryo.TL = Tail.Length..cm..x,
         Embryo.Sex = Sex.x,
         Female.SVL = Snout.to.Vent..cm..y,
         Female.TL = Tail.Length..cm..y,
         Mother.Wt = Weight..g.
         )%>%
  # Removing erroneous data point
  filter(Female.SVL > 50)

Increase in clutch size with female age

re_clutch <- re%>%
  select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
  group_by(Field.Code)%>%
  summarise(clutch.size = n(),
            Clutch.wt = sum(Egg.Weigth..g..),
            Mother.wt = last(Mother.Wt),
            Female.SVL = last(Female.SVL))%>%
  mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))

re_clutch%>%
  ggplot(aes(Female.SVL, clutch.size))+
  geom_point(size =3)+
  geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
  labs(x = "Female SVL (cm)", y = "Clutch Size")

Generalised linear model to test change in clutch size with female SVL:

Clutch size is poisson distributed with mean 8.1052632 and variance 10.3216374

mcfadden <- function(x){
  
  r2 <- 1 - (x$deviance/x$null.deviance)
  
}

re_clutch%>%
  select(Female.SVL, clutch.size)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance),
         r.squared =  map_dbl(stat, ~mcfadden(.)))%>%
  select(sumr, stat,r.squared)%>%
  unnest()%>%
  select(c(term:p.value, r.squared))
term estimate std.error statistic p.value r.squared
(Intercept) -0.3764929 0.7174214 -0.5247862 0.5997318 0.560094
Female.SVL 0.0256926 0.0073002 3.5194451 0.0004325 0.560094

The number of eggs borne by females increase by 1.0253151 for every 1 cm increase in female SVL.

Change in overall reproductive effort with age

re_clutch%>%
  ggplot(aes(Female.SVL, rcm))+
  geom_point(aes(col = clutch.size), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed")+
  scale_y_continuous(name = "Relative clutch mass")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Clutch size")+
  theme(legend.text = element_text(size = 10))

Linear model to determine relationship between relative clutch mass and female SVL:

re_clutch%>%
  select(Female.SVL, rcm)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(rcm ~ Female.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.6886163 0.4077131 1.6889727 0.1170168 0.063695
Female.SVL -0.0038362 0.0042458 -0.9035133 0.3840327 0.063695

The overall reproductive effort does not change with female age.

Change in reproductive effort per embryo with female age

Does the effort per embryo change with female age?

# calculating reproductive effort per embruo

re_embryo <- re%>%
  select(Field.Code, Embryo.Code, Egg.Weigth..g..,  Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
  group_by(Field.Code)%>%
  mutate(clutch.wt = sum(Egg.Weigth..g..), 
         Female.wt = Mother.Wt - clutch.wt)%>%
  group_by(Field.Code, Embryo.Code)%>%
  summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
            Female.SVL = last(Female.SVL),
            Embryo.Sex = Embryo.Sex,
            Embryo.SVL = Embryo.SVL)%>%
  ungroup()

# using residuals to control for embryo development

emsvlinv <- lm(inv ~ Embryo.SVL, data = re_embryo)

# Plotting

re_embryo%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")

Linear model to detemine relationship between female SVL and investment per embryo:

We have controlled for the effect of embryo development.

re_embryo%>%
  select(Female.SVL, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Female.SVL + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.0858802 0.0260524 3.296444 0.0017188 0.4433341
Female.SVL -0.0007894 0.0002588 -3.050593 0.0035109 0.4433341
Embryo.SVL 0.0023248 0.0004171 5.573820 0.0000008 0.4433341

The relative egg mass (controlled for embryo development) reduces with female SVL.

Difference in reproductive effort for male and female embryos

Does female investment differ by sex of the embryo?

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  modelr::add_residuals(emsvlinv)%>%
  ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
  geom_point(aes(col = Embryo.SVL), size = 3)+
  geom_smooth(method = "lm", linetype = "dashed", size = 1)+
  scale_y_continuous(name = "Relative egg mass (residuals)")+
  scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
  scale_color_viridis(name = "Embryo SVL (cm)")+
  scale_shape_discrete(name = "Embryo Sex")

Linear model to test for the difference in female investment by embryo sex:

We controlled for the effect of embryo development.

re_embryo%>%
  filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
  select(Embryo.Sex, inv, Embryo.SVL)%>%
  nest()%>%
  mutate(mod = map(data, ~lm(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
         sumr = map(mod, broom::tidy),
         stat = map(mod, broom::glance))%>%
  select(sumr, stat)%>%
  unnest()%>%
  select(term:p.value, r.squared)
term estimate std.error statistic p.value r.squared
(Intercept) 0.0165286 0.0090610 1.824151 0.0735636 0.3859103
Embryo.SexMale -0.0083114 0.0045803 -1.814613 0.0750384 0.3859103
Embryo.SVL 0.0022764 0.0004455 5.109887 0.0000042 0.3859103

Female investment does not differ significantly with embryo sex.